Simulation of a Composite with a Polyhydroxybutyrate (PHB) Matrix Reinforced with Cylindrical Inclusions: Prediction of Mechanical Properties

Biocomposite development, as a sustainable alternative to fossil-derived materials with diverse industrial applications, requires expediting the design process and reducing production costs. Simulation methods offer a solution to these challenges. The main aspects to consider in simulating composite materials successfully include accurately representing microstructure geometry, carefully selecting mesh elements, establishing appropriate boundary conditions representing system forces, utilizing an efficient numerical method to accelerate simulations, and incorporating statistical tools like experimental designs and re-regression models. This study proposes a comprehensive methodology encompassing these aspects. We present the simulation using a numerical homogenization technique based on FEM to analyze the mechanical behavior of a composite material of a polyhydroxybutyrate (PHB) biodegradable matrix reinforced with cylindrical inclusions of flax and kenab. Here, the representative volume element (RVE) considered the geometry, and the numerical homogenization method (NHM) calculated the macro-mechanical behavior of composites. The results were validated using the asymptotic homogenization method (AHM) and experimental data, with error estimations of 0.0019% and 7%, respectively. This model is valuable for predicting longitudinal and transverse elastic moduli, shear modulus, and Poisson’s coefficient, emphasizing its significance in composite materials research.

Biocomposites offer various advantages due to the renewable origin of their components.For example, they help to reduce dependence on petroleum-based products, promote the use of renewable raw materials, and contribute to the recovery of by-products and agricultural waste [2,[14][15][16][17][18]. Biocomposites are a safe and environmentally sustainable alternative to fossil-derived materials [19,20].Moreover, the biodegradability of these Polymers 2023, 15, 4727 2 of 15 materials aligns with multiple regulations concerning plastic usage and the shift towards a circular economy across various regions, including the Americas, the Caribbean, the European Community, and the United Nations [21][22][23][24].Due to their versatility, biocomposite materials emerge in diverse industry applications, for example, vascular prostheses [25,26], laminated materials [27], and packaging materials [28,29], among other products.
Because of their structural diversity and biodegradability, polyhydroxyalkanoates (PHAs) are exceptional for replacing plastics.PHAs are polymers produced by Grampositive and Gram-negative bacteria from at least 75 genera [30].The most known member of the PHA family is polyhydroxybutyrate (PHB), a highly crystalline, naturally occurring biopolymer [31,32].PHB deficiencies in mechanical properties can be overcome by producing composites with fibers.Natural fibers used as fillers in composites improve characteristics such as stiffness, tensile strength, density, biodegradability, and the cost of the final applications [33].Bio-sourced and biodegradable PHB composites have gained notoriety because various fibers, nanofibers, and modified fibers have been incorporated as PHB fillers, widening the range of applications from packaging to 3D printing [34,35].Yet, the development of biocomposites requires speeding up the design process, reducing production costs, and facilitating the characterization of materials.
Materials simulation methods provide an alternative to the whole process and offer additional advantages, such as non-destructive tests, to characterize the material that significantly reduces the costs [36][37][38], optimizes models through parameterized programming [39,40], and offer the possibility of validating the methodology by comparing results with other methods and experimental data or exploring complex scenarios [41,42].In addition, including different parameters in the simulation methodology, such as volumetric fraction, reinforcement geometries, and orientation, reduces errors due to the overestimation of mechanical properties that occur with some classical theoretical models [43,44].With a validated model and a parametric analysis, it is possible to gain insight into the influence of different parameters and their interactions on the mechanical behavior of composites [45].
The inclusion theory, as formulated in Eshelby's work [46], has been widely used to predict composite materials' overall properties.This approach commonly works with a single inclusion within a sparsely distributed framework in two-phase composites [47].While versatile, the model assumes ellipsoidal inclusions, a representation that may not cover all practical scenarios.It is essential to highlight that the Eshelby method is particularly effective in significantly diluted concentration conditions.However, various alternative methods have been devised, including the Mori-Tanaka approach.This method computes the average internal stress in the matrix of a material with uniformly distributed macroscopic inclusions.The research reveals that the average stress in the matrix remains constant throughout the material, regardless of its position within the domain [48].Additionally, this model has been subject to refinements, which redefined and extended the Mori-Tanaka approach to apply to composite materials, accounting for anisotropic and ellipsoidal phases [49].
Another noteworthy model is the self-consistent method, based on the premise that a particle is embedded in an effective medium with unknown properties [50].This method calculates the effective properties by iteratively adjusting the local properties of each phase until reaching self-consistency.Self-consistency occurs when the average strain in each phase, weighted by the volume fraction of that phase, equals the macroscopic strain.In contrast to the Eshelby and Mori-Tanaka methods, the self-consistent method yields implicit equations for the effective properties, often requiring numerical iteration.
The effective properties of composite materials can also be calculated using a homogenization procedure such as asymptotic homogenization (AH).The AH scheme is widely used to predict effective material properties due to its rigorous mathematical background based on perturbation theory.The boundary value problems of partial differential equations have been studied for application in many diverse engineering fields.The AH method is mathematically rigorous, but it requires a complex and lengthy derivation procedure to obtain effective material properties [51].AH assumes that any field quantity, such as the displacement, can be described as an asymptotic expansion, which, once replaced in the governing equations of equilibrium, allows the effective properties of the composite material to be derived.
Another approach is numerical homogenization (NHM), which employs the finite element method (FEM).This numerical method demonstrates high versatility, adapting well to composites with various geometric inclusion configurations, especially those with intricate fiber arrangements [52].NHM enables efficient calculations of effective properties by utilizing a unit cell model with appropriate periodic boundary conditions.Additionally, NMH goes beyond determining mean-field solutions for highly heterogeneous problems.It also includes the assessment of local fluctuations, a crucial consideration in numerous applications [53].The finite element analysis obtains the local solutions and averages the effective properties of the composite [38,[54][55][56], while statistical tools, such as experimental designs, reduce the number of runs (simulations) required to obtain reliable results [36,45].
This work presents the simulation, using a numerical homogenization technique based on FEM, to analyze the mechanical behavior of a composite material of a polyhydroxyalkanoate (PHA) biodegradable matrix reinforced with cylindrical inclusions of flax.This simulation considers three aspects: (1) An appropriate representation of the geometry, microstructure, and periodicity; for this purpose, a representative volume element (RVE) was used [38,56,57]; (2) Mesh element selection, which establishes adequate boundary conditions that represent the forces that act in the system [58][59][60]; and (3) An efficient numerical method to speed up the calculation time of the simulation [37,61].The results were validated by response surface methodology, AHM, and experimental data.The simulations ran under a fractional experimental design.

Experimental Design
A fractional factorial design (DDF) with three factors and tree levels was used (RStudio 2023.06.2 Build 561, Posit Software, PBC, Boston, MA, USA).The factors were chosen from the literature due to their impact on the mechanical behavior of composites: volumetric fraction (V f ), radius, and aspect ratio between the longitude and radius (ρ) of the reinforcements [62][63][64].The factors and levels are summarized in Table 1.The selected matrix material was polyhydroxybutyrate (PHB), and the reinforcement material was flax fibers.

RVE Generation and Reinforcement Randomization
The RVE and the reinforcement sizing were established within the parameterization section using a random angle vector for orientation.These angles were then used to define the location of the inclusions using rotation matrices.One thousand points on the circumferences of the cylindrical reinforcements were located and validated within the RVE limits.The cell with the inclusions was generated(ANSYS APDL 2022 R1 Build 22.1 Canonsburg, PA, USA), as shown in Figure 1.

Meshing
The discretization of the continuous system was carried out ((SOLID 187 element ANSYS APDL 2022 R1 Build 22.1 Canonsburg, PA, USA), a tetrahedron with ten nodes and three degrees of freedom at each node.The number of elements generated depended on the number and length of reinforcements.

Boundary Conditions
Six local problems were analyzed, three for tensile stress and three for shear stress [65][66][67], summarized in Table 2. Longitudinal tension stress conditions (problem 1) were as follows: A known perpendicular load was applied to the  =  face of the RVE.Likewise, the displacement of the cell in the three planes of symmetry,  = 0,  = 0, and  = 0, was restricted to simulate the continuity of the material [55,68].A homogeneous displacement of the nodes was ensured on the faces  =  , and  =  (faces in which no force was applied).The equivalent procedure was performed for the second and third scenarios.
Conditions of shear force problems (problem 4) were as follows: Two known forces were applied, equal and parallel to the  and  xes, on the faces  =  and  =  .The displacement in the plane of symmetry  = 0 was restricted.Likewise, the movement of the face  = 0 in the direction of x was limited, and that of the face  = 0 in the  direction, to generate antisymmetric to the applied forces.Finally, a homogeneous

Meshing
The discretization of the continuous system was carried out ((SOLID 187 element ANSYS APDL 2022 R1 Build 22.1 Canonsburg, PA, USA), a tetrahedron with ten nodes and three degrees of freedom at each node.The number of elements generated depended on the number and length of reinforcements.

Meshing
The discretization of the continuous system was carried out ((SOLID 187 element ANSYS APDL 2022 R1 Build 22.1 Canonsburg, PA, USA), a tetrahedron with ten nodes and three degrees of freedom at each node.The number of elements generated depended on the number and length of reinforcements.

Boundary Conditions
Six local problems were analyzed, three for tensile stress and three for shear stress [65][66][67], summarized in Table 2. Longitudinal tension stress conditions (problem 1) were as follows: A known perpendicular load was applied to the  =  face of the RVE.Likewise, the displacement of the cell in the three planes of symmetry,  = 0,  = 0, and  = 0, was restricted to simulate the continuity of the material [55,68].A homogeneous displacement of the nodes was ensured on the faces  =  , and  =  (faces in which no force was applied).The equivalent procedure was performed for the second and third scenarios.
Conditions of shear force problems (problem 4) were as follows: Two known forces were applied, equal and parallel to the  and  xes, on the faces  =  and  =  .The displacement in the plane of symmetry  = 0 was restricted.Likewise, the movement of the face  = 0 in the direction of x was limited, and that of the face  = 0 in the  direction, to generate antisymmetric to the applied forces.Finally, a homogeneous

Boundary Conditions
Six local problems were analyzed, three for tensile stress and three for shear stress [65][66][67], summarized in Table 2. Longitudinal tension stress conditions (problem 1) were as follows: A known perpendicular load was applied to the x = A x face of the RVE.Likewise, the displacement of the cell in the three planes of symmetry, x = 0, y = 0, and z = 0, was restricted to simulate the continuity of the material [55,68].A homogeneous displacement of the nodes was ensured on the faces x = A y , and z = A z (faces in which no force was applied).The equivalent procedure was performed for the second and third scenarios.
Conditions of shear force problems (problem 4) were as follows: Two known forces were applied, equal and parallel to the x and y axes, on the faces x = A x and y = A y .The displacement in the plane of symmetry z = 0 was restricted.Likewise, the movement of the face y = 0 in the direction of x was limited, and that of the face x = 0 in the y direction, to generate antisymmetric to the applied forces.Finally, a homogeneous displacement of the nodes on the z = A z faces was ensured.The equivalent procedure was performed for the fifth and sixth scenarios.displacement of the nodes on the  =  faces was ensured.The equivalent procedure was performed for the fifth and sixth scenarios.

RVE Analysis and Solution
Hooke's law (1) was used to integrate the element's nodes, material properties (stiffness), and deformation after applying the loads.

RVE Analysis and Solution
Hooke's law (1) was used to integrate the element's nodes, material properties (stiffness), and deformation after applying the loads.

RVE Analysis and Solution
Hooke's law (1) was used to integrate the element's nodes, material properties (stiffness), and deformation after applying the loads.

RVE Analysis and Solution
Hooke's law (1) was used to integrate the element's nodes, material properties (stiffness), and deformation after applying the loads.

RVE Analysis and Solution
Hooke's law (1) was used to integrate the element's nodes, material properties (stiffness), and deformation after applying the loads.

RVE Analysis and Solution
Hooke's law (1) was used to integrate the element's nodes, material properties (stiffness), and deformation after applying the loads.

RVE Analysis and Solution
Hooke's law (1) was used to integrate the element's nodes, material properties (stiffness), and deformation after applying the loads.
where σ is the stress vector, ε is the strain vector, and C is the stiffness matrix.The matrix C was inverted to obtain the matrix S "compliance."This matrix was used to rewrite the deformations as an unknown dependent change, allowing the deformation to be rewritten in terms of the stress and the elements of the matrix S [66].The engineering constants were calculated with Equation (2).
A numerical homogenization, which relies on finite element simulations utilizing subroutines that calculate volume, deformation, and stress for each RVE element, was subsequently used to compute the engineering constants, as described in Table 3 [55].
Table 3. Compliance element matrix and engineering constants for orthotropic materials [66].

Constant Definition
1

Validation, NHM vs. AHM
The simulation outcomes were compared with the results from the asymptotic homogenization method (AHM) [41,55].Figure 3 shows NHM (blue asterisk) vs. AHM (green squares), revealing that the calculated engineering constants align closely.The mean average discrepancy between the two simulations was 0.0019%.The same behavior was observed with the other constants described in Table 3. mean average discrepancy between the two simulations was 0.0019%.The same behavi was observed with the other constants described in Table 3.

Simulation vs. Experimental Results
To compare the simulation results with experimental data, we used the da generated by Yan et al. from a biodegradable matrix and cylindrical inclusions to matc the parameters included in our research (Table 4) [43].The referred study use transversely isotropic fibers; therefore, it was necessary to rotate the tensioner accordin to the fibers' orientation.
Table 5 shows the NHM simulation results for two  and the error estimation Because the  was not specified, two approximate data were used, taking the Yan et a micrographs as a reference.The  dramatically impacted the results, especially th longitudinal module, probably due to the inclination angles used for the fibers' rotation

Simulation vs. Experimental Results
To compare the simulation results with experimental data, we used the data generated by Yan et al. from a biodegradable matrix and cylindrical inclusions to match the parameters included in our research (Table 4) [43].The referred study used transversely isotropic fibers; therefore, it was necessary to rotate the tensioner according to the fibers' orientation.Table 5 shows the NHM simulation results for two ρ and the error estimations.Because the ρ was not specified, two approximate data were used, taking the Yan et al. micrographs as a reference.The ρ dramatically impacted the results, especially the longitudinal module, probably due to the inclination angles used for the fibers' rotation.

Mechanical Behavior of the Simulated Composite
Figure 4a,b showed that an increase in ρ or V f results in an increment in the modulus E 11 .Additionally, it is possible to identify an interaction between these factors.When ρ =1, the mean value of E 11 is 3686 MPa for V f = 8% and 3790 MPa was the mean value for V f = 12%, representing an approximate 2.7% increase.However, for ρ = 3, the value of E 11 is 4150 MPa for V f = 8%, and 4555 MPa for V f = 12%, indicating a 9% increase in the modulus.This suggests that as ρ increases, the impact of V f on reinforcing the elastic modulus of composite becomes more evident.

Mechanical Behavior of the Simulated Composite
Figure 4a,b showed that an increase in  or  results in an increment in the modulus  .Additionally, it is possible to identify an interaction between these factors.When  =1, the mean value of  is 3686 MPa for  = 8% and 3790 MPa was the mean value for  = 12%, representing an approximate 2.7% increase.However, for  = 3, the value of  is 4150 MPa for  = 8%, and 4555 MPa for  = 12%, indicating a 9% increase in the modulus.This suggests that as  increases, the impact of  on reinforcing the elastic modulus of composite becomes more evident.A similar effect is observed in the case of the modulus  .Figure 4c,d, for  = 1, the change in  from 8% to 12% appears to have an insignificant impact on the stiffness of A similar effect is observed in the case of the modulus G 12 .Figure 4c,d, for ρ = 1, the change in V f from 8% to 12% appears to have an insignificant impact on the stiffness of the composite.However, when ρ = 3, the increase in V f results in a 7% stiffness increment.Regarding the coefficient v 12 .In Figure 4e,f, an increase in ρ leads to an approximately 1% decrease in v 12 , which aligns with the increased composite's stiffness.In the case of v 12 , the interaction effect of the factors is less pronounced, but there is greater variability in the data.The variability of the factor can be examined through the box plot presented in Appendix A.
the composite.However, when  = 3, the increase in   results in a 7% stiffness increment.Regarding the coefficient  12 .In Figure 4e,f, an increase in  leads to an approximately 1% decrease in  12 , which aligns with the increased composite's stiffness.In the case of  12 , the interaction effect of the factors is less pronounced, but there is greater variability in the data.The variability of the factor can be examined through the box plot presented in Appendix A.

Linear Model and Contour Plot to Predict Elastic Modulus
The statistical analysis of the simulation results appears in Figure 5.A Pareto diagram (Figure 5a) highlighted the influence of the Young's Modulus factor.Among these factors,  and   were the most significant contributors, and notably, their interaction also had a positive effect on the elastic modulus.In contrast, the Ratio(E_fiber/E_matrix) and the radius size exhibited comparatively less impact on Young's Modulus.The prediction model for Young's Modulus is expressed by Equation ( 3) and excludes the radius factor due to its little impact on the response variable.and ρ due to significant influences.The   -axis indicates that an increment in   had a slight effect on Young's Modulus, but as  increased, the influence of the   became more pronounced.The response surface pointed to the optimal approach for enhancing the elastic modulus.The color map (blue-yellow) indicates the direction of improvement.

Discussion
In this work, we modeled the mechanical behavior of a composite material of a polyhydroxyalkanoate (PHA) biodegradable matrix reinforced with cylindrical inclusions of flax using the NHM.We validated our results with the AHM method for its accurate prediction of macroscopic behaviors [69,70].This characteristic is pivotal, especially when comparing simulation outcomes.AHM has proven helpful in the simulation work of short fiber composites compared with NHM in three-dimensional cell cases [71].
For a second validation, we compared our NHM results with published experimental data, with errors ranging from 0.33% to 6.45% obtained for the transverse and longitudinal moduli, with a noticeable increase in error as  increased.Several aspects could contribute to this discrepancy.First, the simulations assumed perfect contact between the The prediction model for Young's Modulus is expressed by Equation ( 3) and excludes the radius factor due to its little impact on the response variable.where: Young's Modulus (E 11 ) ; Volumetric Fraction V f ; Aspect Ratio (ρ); Modulus Ratio (β).
Figure 5b reveals insightful interactions between the factors.This graph used V f and ρ due to significant influences.The V f -axis indicates that an increment in V f had a slight effect on Young's Modulus, but as ρ increased, the influence of the V f became more pronounced.The response surface pointed to the optimal approach for enhancing the elastic modulus.The color map (blue-yellow) indicates the direction of improvement.

Discussion
In this work, we modeled the mechanical behavior of a composite material of a polyhydroxyalkanoate (PHA) biodegradable matrix reinforced with cylindrical inclusions of flax using the NHM.We validated our results with the AHM method for its accurate prediction of macroscopic behaviors [69,70].This characteristic is pivotal, especially when comparing simulation outcomes.AHM has proven helpful in the simulation work of short fiber composites compared with NHM in three-dimensional cell cases [71].
For a second validation, we compared our NHM results with published experimental data, with errors ranging from 0.33% to 6.45% obtained for the transverse and longitudinal moduli, with a noticeable increase in error as ρ increased.Several aspects could contribute to this discrepancy.First, the simulations assumed perfect contact between the matrix and reinforcement, which may not hold in experimental conditions [72][73][74].Second, it was assumed that the reinforcements had uniform geometry, a simplification that does not align with real-world variations.Finally, the simulated material did not account for aggregation.All these aspects could introduce additional rigidity to the composite material.Future research will enhance the algorithm by incorporating a probability distribution to simulate agglomerations and geometric variations.Nonetheless, the error percentage falls within the range reported by other algorithms for predicting mechanical properties.Some studies suggest that error prediction is due to non-linear elastic behavior or estimation errors when the mesh size is too small.Despite the limitations of the present model, we considered it worthwhile to obtain a good approximation of the compound and to reduce experimental costs [36,71,[75][76][77][78][79].
The significance of the ρ of reinforcement was underscored in our statistical analysis, a finding that aligns with various experimental studies.For instance, in the case of a PLC-hemp mixture, it was observed that Young's Modulus increased by 48% as the ρ of the reinforcements changed from 19 to 26.However, when this ratio exceeded 30, a saturation effect occurred, likely attributable to moisture absorption by the fibers [80].Additionally, another study [81] supports the notion that particles with larger aspect ratios are more effective in enhancing composite stiffness.In their study, a ρ of 10, combined with a volume fraction of 30%, resulted in a 4.5% increase in Young's Modulus and a 1.3% decrease in Poisson's coefficient.Nevertheless, they also observed a saturation effect when the ρ exceeded 10, suggesting the occurrence of a "shear lag" phenomenon that limits stress transmission from the matrix to the fibers.In a PHBV-potato mixture, researchers have noted that fibers with small aspect ratios tend to function more as fillers than reinforcements [82].
When the matrix is reinforced with long particles or fibers and the fiber length is lower than the length of RVE, the matrix and fiber interact as springs connected in series, resulting in a high reinforcement effect.However, if the fiber length exceeds the dimensions of the RVE, the interactions resemble parallel springs, leading to a diminished reinforcement effect.The β and the ρ also interact, and when β is >10, these interactions become dominating factors in determining the directional effective Young's Modulus [83].
Of interest, Robinson et al. [84] highlighted that the effectiveness of fiber length in reinforcing a material depends not only on its diameter but also on the ratio between the elastic modulus of the reinforcement and the resin (β), the Poisson coefficient, and the volume fraction.This relationship is mathematically presented in Equation (4).
where: The fiber diameter (d), the modulus fiber (E f ), the matrix modulus (E m ), the matrix poisson (v m ), the volumetric fraction (V f ).According to this relationship, the minimum length of the reinforcements should be 6 mm; however, the longest inclusion used in the simulation is 2.8 mm, because of which the elastic modulus of the matrix increased between 17 and 31% and the shear by 15-24%.The effect of the β was also analyzed, but for the simulated cases, the effect of the β is not so significant.
Additionally, most studies concur on the influence of ρ on the Elastic Modulus, but there are notable exceptions.In the epoxy-nanocarbon fibers research, the authors reported that reinforcements with a lower ρ led to a greater apparent Young's Modulus [85], attributed to enhanced interfacial adhesion.Also, the authors reported that reinforcements with a lower ρ led to a greater apparent Young's Modulus.This phenomenon results from enhanced interfacial adhesion due to better dispersion and alignment of short fibers.The researchers also emphasized that the V f had a more significant impact than the ρ opposite to the simulation results [85].Furthermore, another study analyzing a PHBV-cellulose composite [86] showed that weak interfacial adhesion softened the matrix modulus, particularly when employing reinforcements with a small ρ.
The simulation works have also analyzed the reinforcement orientation; unidirectional fibers increase longitudinal modulus by 26%, but the transverse modulus decreases to 4.75.This effect depends on interacting with the V f [87].Other work states that a simulation with random fibers is closer to the experimental data [88].The proposed methodology uses random reinforcement orientation to simulate the injected composite material.In the experimental validation, the fiber rotation was limited to −30 • to 30 • since the micrographs of the work used as a reference showed limited fiber rotation [43].The last methodology step presented in this research proposed a non-linear regression model, Equation (3), to predict the Young's Modulus of the compound.The model is considered helpful in predicting the mechanical properties of the composite when parameters such as ρ and V f are varied.For example, increasing ρ to 9.5 and V f to 20% (in coding units) and then substituting into Equation (3) while keeping the β constant estimates the Young's Modulus (E 11 ) to be 7515.2MPa.
It is necessary to consider that the model is a generalization and may incorporate noise.Additionally, the model assumes a constant relationship between the factor and E 11 , which can lead to inaccurate predictions when the factor's level is changed to values significantly distant from those simulated.Non-linear regression models are sensitive to the initial parameter values.
Although this simulation methodology might seem simple, it allows precise predictions of the mechanical behavior of composites, including changes in elastic and shear modulus, in scenarios involving fibers with approximate cylindrical shapes.The regression model enhances result reliability, leading to reduced computational costs, design, and production times for materials with similar characteristics (β,ρ,V f ).For example, it is possible to utilize the model to predict the elastic modulus of a PHB-Jute composite with V f = 16%, β = 7.7, and ρ = 5, resulting in E 11 = 4800 MPa (approximately).Worth noting is that the structure of the composites is complex and requires the design of complex algorithms that include random distributions of the reinforcements in the unit cells, include movement restrictions of the nodes and planes, and thus represent the reality of a composite material.

Conclusions
Many fibers are industrial waste and are underused.Biodegradable matrix composites with fiber reinforcements have been developed because they could potentially replace polymers in olefins.Composites of PHB with natural fibers have great potential for industrial applications.The algorithm developed here facilitated the simulation of the reinforcement of a PHB matrix with cylindrical fiber inclusions and produced reliable results.The algorithm's parameters helped simulate different configurations of the materials.This algorithm relied on the geometric depiction of the material and used numerical homogenization methods for analysis.This model is a valuable tool for predicting both longitudinal and transverse elastic moduli, shear modulus, and Poisson's coefficient.Furthermore, a linear regression model was developed, simplifying the prediction of mechanical properties for alternative material configurations while mitigating the need for extensive computational resources.This simulation methodology sets the basis for future work, including challenging conditions such as porous cells, aggregations, and inclusions with different aspect ratios.This advancement aims to enhance the representation of composite materials in a simulation process.

Problem 2 :
Tensile stress.Axis = y Polymers 2023, 15, x FOR PEER REVIEW 5 of 16displacement of the nodes on the  =  faces was ensured.The equivalent procedure was performed for the fifth and sixth scenarios.

Problem 3 :
Tensile stress.Axis = z Polymers 2023, 15, x FOR PEER REVIEW 5 of 16displacement of the nodes on the  =  faces was ensured.The equivalent procedure was performed for the fifth and sixth scenarios.

Problem 4 :
Shear stress plane = xyPolymers 2023, 15, x FOR PEER REVIEW 5 of 16 displacement of the nodes on the  =  faces was ensured.The equivalent procedure was performed for the fifth and sixth scenarios.

z
= 0. y = 0, restrict movement in direction x.x = A x , y = A y z = A z. x = 0, restrict movement in direction y.Problem 5: Shear stress plane = yz Polymers 2023, 15, x FOR PEER REVIEW 5 of 16displacement of the nodes on the  =  faces was ensured.The equivalent procedure was performed for the fifth and sixth scenarios.

x
= 0. y = 0, restrict movement in direction z. y = A y z = A z .x = A x .z = 0, restrict movement in direction y.Problem 5: Shear stress plane = xz Polymers 2023, 15, x FOR PEER REVIEW 5 of 16displacement of the nodes on the  =  faces was ensured.The equivalent procedure was performed for the fifth and sixth scenarios.

Figure 3 .
Figure 3.Comparison between NHM and AHM.In (a,c,e), the x-axis represents the .In (b,d,f), th x-axis represents the percentage of reinforcements included in the matrix,  .

Figure 3 .
Figure 3.Comparison between NHM and AHM.In (a,c,e), the x-axis represents the ρ.In (b,d,f), the x-axis represents the percentage of reinforcements included in the matrix, V f .

Figure 4 .
Figure 4. Comparison between NHM and AHM.In (a,c,e), the x-axis represents the , and the color bar represents  .In (b,d,f), the x-axis represents the percentage of reinforcements included in the matrix  and the color bar represents .

Figure 4 .
Figure 4. Comparison between NHM and AHM.In (a,c,e), the x-axis represents the ρ, and the color bar represents V f .In (b,d,f), the x-axis represents the percentage of reinforcements included in the matrix V f and the color bar represents ρ.

Figure 5 .
Figure 5. (a)  11 Pareto plot showing the positive impact of the analyzed factors on elastic modulus; (b)  11 response surface.

Figure 5 .
Figure 5. (a) E 11 Pareto plot showing the positive impact of the analyzed factors on elastic modulus; (b) E 11 response surface.

Table 2 .
Boundary conditions for tensile and shear stress.

Table 2 .
Boundary conditions for tensile and shear stress.

Table 2 .
Boundary conditions for tensile and shear stress.

Table 2 .
Boundary conditions for tensile and shear stress.

Table 2 .
Boundary conditions for tensile and shear stress.

Table 2 .
Boundary conditions for tensile and shear stress.

Table 2 .
Boundary conditions for tensile and shear stress.

Table 4 .
Experimental data were used for simulation.

Table 4 .
Experimental data were used for simulation.

Table 5 .
Experimental vs. simulation results by NHM.

Table 5 .
Experimental vs. simulation results by NHM.